///@file Tensor4.cpp
///@brief Implementation of the class representing a fourth order tensor
///@author Arnaud Duval
///@version 0.0
///@date 2010/02/26


#include "Tensor4.h"

///@brief Default constructor
Tensor4::Tensor4()
{
	size = 0;
}

///@brief Constructor giving the dimension
///@param s size of the tensor
Tensor4::Tensor4(unsigned int s)
{
	assert(s > 0);
	size = s;
	t = new Tensor3[s];
	for(unsigned int i = 0 ; i < s ; i++)
	{
		t[i] = Tensor3(s);
	}	
}

///@brief Copy constructor
///@param tens Tensor4 object to copy
Tensor4::Tensor4(const Tensor4 &tens)
{
	size = tens.size;
	t = new Tensor3[size];
	for(unsigned int i = 0 ; i < size ; i++)
	{
		t[i] =tens.t[i];
	}
}

///@brief Copy constructor with temporary object
///@param tens Temporary Tensor4 object to copy
Tensor4::Tensor4(Tensor4 *tens)
{
	size = (*tens).size;
	t = new Tensor3[size];
	for(unsigned int i = 0 ; i < size ; i++)
	{
		t[i] = (*tens).t[i];
	}
}

///@brief Destructor
Tensor4::~Tensor4()
{
	if(size)
		delete[] t;
	size = 0;
}

///@brief Return an identity
///@param dim Dimension of the tensor to return (only dimension 3 case is implemented)
Tensor4 Tensor4::Identity(unsigned int dim)
{
    assert(dim == 3);
    Tensor4 res(3u);
    Tensor2 delta;
    delta = delta.Delta(3u);
    for(int i = 1 ; i <= 3 ; i++)
        for(int j = 1 ; j <= 3 ; j++)
            for(int k = 1 ; k <= 3 ; k++)
                for(int l = 1 ; l <= 3  ; l++)
                    {
                        res(i,j,k,l) = 0.5 * (delta(i,k)*delta(j,l)  +  delta(i,l)*delta(j,k));
                    }

    return res;
}

///@brief Overload of () operator
double& Tensor4::operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
{
	assert((i > 0)&&(i <= size));
	assert((j > 0)&&(j <= size));
	assert((k > 0)&&(k <= size));
	assert((l > 0)&&(l <= size));
	return (t[i-1])(j,k,l);
}

///@brief Overload of = operator
///@param tens Tensor4 object to copy
Tensor4& Tensor4::operator = (const Tensor4 tens)
{
	assert (tens.size != 0);
	size = tens.size;
	t = new Tensor3[size];
	for(unsigned int i = 0 ; i < size ; i++)
		t[i] = tens.t[i];
	return (*this);
}


///@brief Overload of = operator for temporary objects
///@param tens Pointer to Tensor4 object to copy
Tensor4& Tensor4::operator = (Tensor4* tens)
{
	assert ((*tens).size != 0);
	size = (*tens).size;
	t = new Tensor3[size];
	for(unsigned int i = 0 ; i < size ; i++)
		t[i] = (*tens).t[i];
	return (*this);
}

///@brief Overload of += operator
///@param tens Third order tensor to add to current one
void Tensor4::operator += (const Tensor4& tens)
{
	assert(size == tens.size);
	for(unsigned int i = 0 ; i < size ; i++)
		t[i] += tens.t[i];
}


///@brief Overload of -= operator
///@param tens Third order tensor to substract to current one
void Tensor4::operator -= (const Tensor4& tens)
{
	assert(size == tens.size);
	for(unsigned int i = 0 ; i < size ; i++)
		t[i] -= tens.t[i];
}

///@brief Overload of || operator : twice contracted product : A_ijkl B_kl
///@param t Tensor2 right hand operand
Tensor2 Tensor4::operator||(Tensor2 tens)
{
    assert(size == tens.Dim());
    Tensor2 res(size);
    for(unsigned int i = 1 ; i <= size ; i++)
        for(unsigned int j = 1 ; j <= size ; j++)
        {
            res(i,j) = 0.;
            for(unsigned int k = 1 ; k <= size ; k++)
                for(unsigned int l = 1 ; l <= size ; l++)
                    res(i,j) += t[i-1](j,k,l)*tens(k,l);
        }
    return res;
}

///@brief Overload of || operator : twice contracted product : A_ijkl B_kl
///@param A left operand : Tensor4
///@param B right operand : Tensor2
Tensor2 operator ||(const Tensor4& A, Tensor2& B)
{
	assert(A.size == B.Dim());
 	Tensor2 res(A.size);
	for(unsigned int i = 1 ; i <= A.size ; i++)
		for(unsigned int j = 1 ; j <= A.size ; j++)
		{
			res(i,j) = 0.;
			for(unsigned int k = 1 ; k <= A.size ; k++)
				for(unsigned int l = 1 ; l <= A.size ; l++)
					res(i,j) += A.t[i-1](j,k,l)*B(k,l);
		}	
	return res;
}

///@brief Overload of || operator : twice contracted product : A_ij B_ijkl
///@param A left hand operand : Tensor2
///@param B right hand operand : Tensor4
Tensor2 operator ||(Tensor2& A, const Tensor4& B)
{
	assert(A.Dim() == B.size);
	Tensor2 res(B.size);
	for(unsigned int k = 1 ; k <= B.size ; k++)
		for(unsigned int l = 1 ; l <= B.size ; l++)
		{
			res(k,l) = 0.;
			for(unsigned int i = 1 ; i <= B.size ; i++)
				for(unsigned int j = 1 ; j <= B.size ; j++)
					res(k,l) += A(i,j) * B.t[i](j,k,l);
		}
	return res;
}

///@brief Overload of || operator : twice contracted product : A_ijkl B_kl
///@param A Left hand operand : Tensor4
///@param B Right hand operand : SymTensor2
///@warning The fact that the SymTensor is symmetric is not taken into account : should be improved
Tensor2 operator ||(const Tensor4& A, SymTensor2& B)
{
	assert(A.size == B.Dim());
 	Tensor2 res(A.size);
	for(unsigned int i = 1 ; i <= A.size ; i++)
		for(unsigned int j = 1 ; j <= A.size ; j++)
		{
			res(i,j) = 0.;
			for(unsigned int k = 1 ; k <= A.size ; k++)
				for(unsigned int l = 1 ; l <= A.size ; l++)
					res(i,j) += A.t[i-1](j,k,l)*B(k,l);
		}	
	return res;
}

///@brief Overload of || operator : twice contracted product : A_ij B_ijkl
///@param A Left hand operand : SymTensor2
///@param B Right hand operand : Tensor4
///@warning The fact that the SymTensor is symmetric is not taken into account : should be improved
Tensor2 operator ||(SymTensor2& A, const Tensor4& B)
{
	assert(A.Dim() == B.size);
	Tensor2 res(B.size);
	for(unsigned int k = 1 ; k <= B.size ; k++)
		for(unsigned int l = 1 ; l <= B.size ; l++)
		{
			res(k,l) = 0.;
			for(unsigned int i = 1 ; i <= B.size ; i++)
				for(unsigned int j = 1 ; j <= B.size ; j++)
					res(k,l) += A(i,j) * B.t[i](j,k,l);
		}
	return res;
}
///@brief Overload of * operator : product A_ij B_kl
///@param A left hand operand : Tensor2
///@param B right hand operand : Tensor2
Tensor4 operator *(Tensor2& A, Tensor2& B)
{
	assert(A.Dim() == B.Dim());
	Tensor4 res(A.Dim());
	for(unsigned int i = 1 ; i <= A.Dim() ; i++)
		for(unsigned int j = 1 ; j <= A.Dim() ; j++)
			for(unsigned int k = 1 ; k <= A.Dim() ; k++)
				for(unsigned int l = 1 ; l <= A.Dim() ; l++)
					res.t[i-1](j,k,l) = A(i,j) * B(k,l);
	return res;
}

///@brief Overload of + operator : addition of 2 Tensor4
///@param t1 left hand operand : Tensor4
///@param t2 right hand operant : Tensor4
Tensor4 operator +(const Tensor4& t1, const Tensor4& t2)
{
	assert(t1.size == t2.size);
	Tensor4 res(t1.size);
	for(unsigned int i = 0 ; i < t1.size ; i++)
		res.t[i] = t1.t[i] + t2.t[i];
	return res;
}

///@brief Overload of - operator : substraction of 2 Tensor4
///@param t1 left hand operand : Tensor4
///@param t2 right hand operand : Tensor4
Tensor4 operator -(const Tensor4& t1, const Tensor4& t2)
{
	assert(t1.size == t2.size);
	Tensor4 res(t1.size);
	for(unsigned int i = 0 ; i < t1.size ; i++)
		res.t[i] = t1.t[i] - t2.t[i];
	return res;
}

///@brief Overload of * operator : product Tensor4 - scalar
///@param t left hand operand : Tensor4
///@param s right hand operand : scalar
Tensor4 operator *(const Tensor4& t, const double& s)
{
	Tensor4 res(t.size);
	for(unsigned int i = 0 ; i < t.size ; i++)
		res.t[i] = t.t[i] * s;
	return res;
}

///@brief Overload of * operator : product scalar - Tensor4
///@param s left hand operand : scalar
///@param t right hand operand : Tensor4
Tensor4 operator *(const double& s, const Tensor4& t)
{
	Tensor4 res(t.size);
	for(unsigned int i = 0 ; i < t.size ; i++)
		res.t[i] = t.t[i] * s;
	return res;
}

///@brief Overload of / operator : division by a scalar
///@param t left hand operand : Tensor4
///@param s right hand operand : scalar
Tensor4 operator /(const Tensor4& t, const double& s)
{
	assert(s != 0.);
	Tensor4 res(t.size);
	for(unsigned int i = 0 ; i < t.size ; i++)
		res.t[i] = t.t[i] / s;
	return res;
}

///@brief A function to compute cross product A_ij B_kl because of a problem with operator *
///@param A First Tensor2 A_ij
///@param B second Tensor2 B_kl
Tensor4 Tensor4::CrossProduct(Tensor2 A, Tensor2 B)
{
    assert(A.Dim() == B.Dim());
    Tensor4 res(A.Dim());
    for(unsigned int i = 1 ; i <= A.Dim() ; i++)
        for(unsigned int j = 1 ; j <= A.Dim() ; j++)
            for(unsigned int k = 1 ; k <= A.Dim() ; k++)
                for(unsigned int l = 1 ; l <= A.Dim() ; l++)
                    res.t[i-1](j,k,l) = A(i,j) * B(k,l);
    return res;
}


///@brief Overload of << operator
///@param s out stream
///@param t Tensor4 object to put in the out stream
std::ostream& operator<<(std::ostream &s, const Tensor4& t)
{
	for(unsigned int i = 0 ; i < t.size ; i++)
	{
		s << t.t[i] << "\n";
		if(i != t.size -1)
			s << "\n";
	}
	return s;
}
